### RD functions ###

#### * FUNCTION * rd coefficients ####
rdmodel <- function(datause,
                    depvar,
                    bw=10^10,
                    cont="",
                    treat="until1789",
                    dist="distestate",
                    clust="generalite",
                    loc="nearET + latitude + longitude",
                    summary=T){
  
  datause <- as.data.table(datause)
  
  dta <- datause[apply(datause[,..dist],1,min)<bw,]
  f <- as.formula(paste("scale(",depvar,") ~ ", treat,"+",loc,"+",paste(cont,collapse="+"),"| 0 | 0 | generalite",sep=""))
  model <- felm(f, data=dta)
  #if(summary==T){
  out <- model %>% summary %>% .$coefficients %>% .[treat,]
  #}
  if(summary==F){
    out <- model 
  }
  return(out)
}


rdplot <- function(dat,var,labvar=var,exclude="",cont=controls,main="",cexaxis=1.3, ylim=c(-1,1)){

  coefs <- lapply(xs*1000, FUN=rdmodel, datause=dat, cont=cont, loc="latitude + longitude + latitude*longitude", depvar=var)
  coefs <- do.call(rbind,coefs)
  plot(xs,coefs[,1], ylim=ylim, ylab="",yaxt='n',main=main, cex.main=1.5, xlab="")
  axis(4,at=c(-1,0,1))
  axis(2,at=0,label=labvar,las=1,tick=F,cex.axis=cexaxis)
  segments(xs,coefs[,1]+1.96*coefs[,2],xs,coefs[,1]-1.96*coefs[,2])
  abline(h=0,lty=2)
  
  
}


getmodel<-function(datause,outcome,treat,clust,controls,scale=T){
  
  if(controls[1]!=""){
    f <- as.formula(paste("depvar ~ ",treat,"+",paste(controls,collapse="+")))
  }
  if(controls[1]==""){
    f <- as.formula(paste("depvar ~ ",treat))
  }
  
  depvar<-as.numeric(datause[,outcome])
  if(scale==T){
    depvar <- scale(depvar)
  }
  lm <- lm(f, data=datause)
  
  N<-length(lm$fitted.values)
  if(length(clust)>1){
    vcov<-cluster.vcov(lm, datause[,clust])
    model<-coeftest(lm,vcov)
  }
  if(length(clust)==1 & clust[1]!=""){
    vcov<-cluster.vcov(lm, datause[,clust])
    model<-coeftest(lm,vcov)
  }
  if(clust[1]==""){
    model<-coeftest(lm)
  }
  return(model)
}

